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Introduction 

By expanding the solution to the initial value problem 

I)' = f(t,y) y(t 0 ) = y 0 (l.la,b) 

in a Taylor series about t 0 , one obtains a local solution which is valid within its radius of 
convergence R 0 . If the series is evaluated at < 1? where t l < R 0 , one obtains an approx- 
imation to y(ti), and the solution may then be expanded in a new series about t x . The 
solution may of course then be extended to the point and so forth, so that by a process 
of “analytic continuation” one obtains a piecewise polynomial solution to (1.1). 

When the derivatives of y can be obtained cheaply, the method of Taylor series, or 
analytic continuation, offers several advantages over other methods. Foremost among these 
is that a Taylor series solution provides information that is not provided by other methods. 
This includes derivative information, local radius of convergence, and the location of poles 
in the complex plane. Another advantage is that both the stepsize and order can be easily 
changed, so that an optimal stepsize can always be chosen. Finally, Taylor series integration 
provides a piecewise continuous solution to the ODE, so that it is never necessary to 
interpolate at intermediate points. 

Since analytic continuation is efficient only when high order derivatives can be ob- 
tained inexpensively, much effort has been made to develop recursive techniques that can 
generate the derivatives numerically. The first known attempt to incorporate recursion 
into the Taylor series method was by J.R. Airey in 1932 [1]. Miller was the first to employ 
a full recurrence scheme in which all series coefficients were calculated numerically, and 
the first to develop a systematic approach which was applicable to a class of differential 
equations [2,3]. Other early work includes that of Gibbons [4] and Moore [5], who devel- 
oped programs that automatically generate the series coefficients for differential systems 
in canonical form. Other work concerning automated differentiation is that of Leavitt [6], 
Nikolaev [7], Barton et. al [8,9], Norman [10], Kedem [11], Rail [12], Chang [13], Corliss 
and Chang [14], Irvine and Savageau [15], and Savageau and Voit [16]. 

Of the above, the most complete and systematic treatments of the Taylor series method 
are those of Barton, Willers, and Zahar [8,9], Chang [13,17,18] and Corliss and Chang [14], 
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and Irvine and Savageau [15]. Barton et.al. were the first to develop and implement a 
general purpose algorithm that automatically reduces a given differential system to canon- 
ical form and generates the series coefficients without any previous manipulations by the 
user. Chang and Corliss developed a similar approach along with computer programs that 
are fully portable. Their programs also provide information on the location and order of 
primary singularities. Irvine and Savageau developed highly efficient series methods for 
differential equations that can be expressed in the form of an “S-system” . Their approach 
obviates the need to generate problem-dependent, source code for each new" problem. 

Corliss and Chang and Irvine and Savageau compared their codes with several stan- 
dard ODE solvers on the problem set of Hull et.al. [19]. Compared to the IMSL routines 
DVERK (Runge-Kutta method) and DGEAR (variable order Adams method), Corliss 
and Chang report that their code was significantly faster for error tolerances in the range 
10~ 9 to 10 -12 . They conclude that the Taylor series method is superior to the variable 
order Adams method and the Runge-Kutta method when stringent accuracies are required 
and when solving small systems. Irvine and Savageau compared their method with two 
standard Runge-Kutta schemes, a variable order Adams method, and Gear’s method. On 
an S-system from cellular biology, they report solution times that were significantly faster 
than all four methods for error tolerances ranging from 10 -1 to 10 -15 . In some cases, 
their code was faster by one or two orders of magnitude. On Hull’s problem set, their 
results compared favorably both in terms of solution times and number of function eval- 
uations. For small error tolerances, their code was up to 20 times faster. One concludes 
that the Taylor series method is indeed competitive with other methods, especially when 
high accuracies are required. 

In this paper we propose an extension of the Taylor series method which improves its 
accuracy and stability while also increasing its range of applicability. The basic idea is 
to formulate the method as a generalized implicit method. Specifically, given an interval 
[f n ,f n+ i], we consider a series solution to (1.1) which is expanded about the intermediate 
point t n + /i h, where h is the stepsize and // is an arbitrary parameter called an expansion 
coefficient. 

The main results of the paper are the following, (i) We derive the discretization error 
for a series solution to (1.1) which is expanded about the point t n + yh. It is shown that 
a series of degree k has exactly k expansion coefficients which raise its order of accuracy. 
The increase is one order if k is odd, and two orders if k is even. For k odd, it will be seen 
that a class of Turan type, one-stage implicit Runge-Kutta methods are recovered [20]. (ii) 
We consider stability for the problem y' = Xy, Re (A) < 0, and identify several A-stable 
schemes, (iii) We show that, for all series of degree three or higher, local extrapolation 
can be used to increase accuracy by two additional orders, (iv) We discuss variable step 
schemes, (v) We present numerical results for a range of problems, including ODE’s with 
a singular point and stiff systems. Our results indicate that implicit Taylor series methods 
are an effective integration tool for most problems. 
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Taylor Series of Degree One 


We assume that f(t, y ) is a real-valued function which is analytic at the point ( t n , y n ), 
and that the problem 

y — fif-y), y{tn) = yni (2.ia,b) 

has a unique, continuous solution which can be expressed as a power series that converges 
for all t in the interval \t — t n \ < R n . We denote the stepsize / n +i — tn by h, where 
h < R n . 

A Taylor series of first degree expanded about t n gives the equation 

h 2 

y{tn+ 1) = yitn ) T hf{t n ,yn) + f (£rn y{^n))i t n < t n + 1, (2-2) 

from which one obtains Euler’s method 


Vn+ 1 — Vn "h hf(t, n , y n ) (2-3) 

and its associated discretization error tL. f f (£ n , y(£, n )) y''{£,n)- The error may also be 

expressed by the expansion 

e n+ 1 := y(t n + 1 ) - y n + 1 = ~2 In + "g"/n + 7^ fn + (2.4) 

where /' := f f (t n ,y n ), and similarly for /" and 

If the series is expanded about t n + 1 , one obtains the backward Euler method 

Vn + 1 == Vn T tl /(tn+1 , ?/n+l ) ? (2.o) 


where ?/„ +1 must be obtained implicitly. By expanding /(t n +i 5 2/n+i) in a Taylor series 
about ( t n ,y n ), one derives from (2.5) the discretization error 


^n+l ■ — y{tn+ 1) ?/n+l 


r h 2 , h 3 

TT fn + ~7T /n + 


L 2 


6 




+ < 3 /»' + 8 ^ + 12 ^ /n>] + °(fc 5 ). (2.6) 

Both the forward and backward Euler methods provide a linear approximation to y 
on I n := [t n ,t n+ 1 ]. This approximation can be written 


yn — yn+fi T (t ^n+^) /(^n+pi; 2/n+^i)i t n < t < ^n+l> (2-7) 


where is zero or one, respectively. 
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If we define t, n+fl := t n + yh, where 0 < y < 1, Eq. (2.7) becomes a generalized 
approximation to y on I n . In this case the approximation will pass through the intermediate 
value y n +fi , which is unknown. By observing that y n must equal y n when t = t n , one can 

see that y n + tl must be determined from 

yn+fi = yn "k y h 2/n+//)* (2.8a) 

Note that this is just the backward Euler method for the partial step from t n to t n+ll . 
Once y n +n is determined from (2.8a), one sets t = £ n +i, y n = y n +i in (2.7) to obtain 

Vn+i = y n +fi + (1 - fi)hf(t n+lt ,y n+tl ). (2.8b) 

One observes that this is just, the forward Euler method for the partial step from t n+/1 to 
£ n+ 1 - Thus, two steps are required to advance the solution to t n+ 1 - a backward Euler 
intermediate step followed by a forward Euler final step. 

To determine an optimal expansion point t n+fl , one obtains the truncation error for 
(2.8). Since (2.8a) is the backward Euler method with stepsize yh, one uses (2.6) to obtain 
an expansion for y n+At , which may then be substituted into (2.8b). One gets from (2.6), 

Un+n = y n + y h In + (yh) 2 f' n + (fn + fn) 

+ (C + 2 S'i + 3 ^ l’„) + 0(A 5 ). (2.9) 

Upon expanding /(f n+/i , y n + fl ) about (£„,?/„), and evaluating powers of (y n+/i - y n ) by 
way of (2.9), one obtains from (2.8b) 


Vn + 1 = yn + hf n + fi h 2 /' + y 2 y (/" + ~ /' ) 

+y’-j tc + 2^/: + 3^/:) + o^. 


( 2 . 10 ) 


The error expansion for (2.8) is then 


h 2 h 3 r r) f 

e n +i := y(t n + 1 ) - Vn+ i = (l-2y)—f' n + — [(l-3// 2 )/" - 3 y 2 /' 


+ £l(i - in 3 )c - y(s + 12 /„)] + o(h s ). 


24 L' 


dy 


( 2 . 11 ) 
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Taking /i — the marching scheme becomes 

2/n-f^- Vn ~2 ) 

h . s 

2/n+l 2/rc+^ ^ /(^n+ ^ 5 ) 5 


with discretization error 


®n+l 


h 3 

= — if" 

24 


Q 9/n f / 

3 % /n) 


(2.12a) 

(2.12b) 


, / fin _ p r" _ 

48 9?/ Jn 


o dfn fl \ 

6 dy }n ’ 


+ 0(h 5 ). 


(2.13) 


The choice y = \ thus leads to a second-order scheme. In view of the central expansion 
point t n+ 1, we will refer to (2.12) as the “centered Euler method”. 

It should be noted that, corresponding to (2.12), the linear expansion (2.7) approx- 
imates y(t ) with an error at most equal to |/ T ' ( | + 0(h 3 ). This is one-fourtli the error 
bound on (2.7) corresponding to the forward and backward Euler methods. 

We should also mention that the explicit halfstep (2.12b) can be equivalently expressed 
as 

y n + 1 = Vn + hf(t n+ i,y n+ i). (2.14) 

When written in this form, it is clear that the centered Euler method is similar to the 
trapezoidal rule 

2/n+l = 2/n "b T [/(^T!) 2/n) T f (tn+l: 2/n+l)] • (2.15) 

It is also similar to the one-stage Gauss implicit scheme [21], [22], 


2/n-|-l — 2/n + ^1 


(2.16a) 


h = h f(t n+ i,y n + Ut), 

which is equivalent to the implicit midpoint rule [22], 


(2.16b) 


Vn+i = Vn + hf(t n+ L,±(y n + y„ + i)). 


(2.17) 


It can be shown in fact that the centered Euler method and the implicit midpoint rule are 
equivalent schemes. One can also show that the trapezoidal linear approximation to y on 
I n : 

, u x \ fi^ni 2/n) 4" f (tn+1 > 2/n+l) (o 18) 

2 MO = 2 Jn + ( t~t n ) , (^-loj 

approximates y with a maximum error equal to that of the centered Euler method, |/(,| + 
0(h 3 ). One concludes that these three schemes offer nearly identical benefits over the 
forward and backward Euler methods. 
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Taylor Series of Degree Two 


Analogous to Eq. (2.7), a Taylor series of second degree expanded about the point 
ln+ii approximates y on I n by 


Vn — l/n+fi T (i tn+n) f (tn+fi'! Un+fj.) T 


(t 


^n+/i) 2 


/ {tn+fii Vn+n)- (3-1) 


As in the previous section, the expansion coefficient y defines a family of marching schemes. 
Choosing //. = 0 and setting t = t n+ i , y n = y n +i, one obtains the Taylor series method of 


order two 


with local error 


Vn + 1 


^n-f-1 


y n + hf(t n ,y n ) + — f'(t n .y n ) 


h 3 h 4 h 5 

n fit _J_ 'J_ fW , n £fW 
6 


f + — f + — f 


+ 


0! 


(3.2) 


(3.3) 


Similarly, choosing y = l and setting t = t n , y n = y n , one obtains the implicit method 


Vn+1 = Vn + hf(t n+ i,y n+ i) — f'(t Tl+ i,iy n+ i). (3.4) 

By expanding both f(t n+ i,y n+ i) and f'(t n+1 ,y n+ i) in a Taylor series about ( t n ,y n ), one 
derives from (3.4) the error expansion 

h 3 h 4 df 

e„ +1 = T /" + fj(3 /”' + 4^0 

+ ^ (6 C + 15 ^ C + 10 ^ O + 0(5 6 ). (3.5) 

When 0 < ju < 1, (3.1) leads to a two-halfstep, implicit/explicit marching scheme, 
analogous to (2.8). Setting t = t n , y n = y n , one obtains the implicit halfstep 

( L \ 2 

Vn+fi = Vn 4 " fi h f Vn+fi) ^ f {tn+fit Vn+fi) • (3.6a) 

Setting t = t n+ 1 , y n = y n +i- one obtains the explicit halfstep 


Vn+ 1 — Vn+fi T" (1 fi) h f , y n +/i) 4~ 


[(1 - fi) hf 


f'{tn+tl,yn+n)- (3.6b) 


An error expansion for (3.6) can be derived by proceeding as in the previous section. 
Utilizing the fact that (3.6a) is just (3.4) for the partial step y h, one obtains an expansion 
for t/n+jj by way of Eq. (3.5). One gets 


Vn+» = Vn + fihfn + 



(p ( fill , cy dfn fll\ 

12 ^ /n + dy Jn) 
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W" + 3 ^»' + 2 ^U») + 


(3.7) 


An expansion for y n+ i is then obtained from (3.6b) by expanding f (t n+ll ,y n+fl ) and 
f'(t n+ll ,y n+ M ) about (t n ,y n ), and evaluating powers of (t/„+ m - Vn) by way of (3.7). Re- 
taining terms through 0(h 5 ), one finds that 


2 /n+l 


Vn + */n + y/; + 3//(l-//)^/" 


2 h 4 

» 24 


( 6 - 8 ,.)/"' - 4 „^/"] 


(3.8) 


+ ,. 3 ^[(l0-15^)/r - 15^^/"' - 10 f/"] + 0(h 6 ). 


Subtracting the left hand side of (3.8) from y(t n+ 1 ) and the right hand side from the Taylor 
series expansion of y(t n + 1 ), one obtains 

5n+i = (1 - 3 y + 3 ,. 2 ) /" + ^ (l - 6 y? + 8 ,. 3 )/"' + 4 ,x 3 /" 


h 5 r 
+ 5f 


24 L 

df n 


(1 - 10^ 3 + 15/u 4 )/"" + 15/^-/"' + 10 y 3 ^fn I + 0(/i 6 ). 


dy 


dy 


Setting the leading coefficient to zero, 

1 — 3 ,. + 3 y? = 0, 


(3.9) 


(3,10) 


one finds 




1 , . \/3 

— i % , 

2 6 


(3.11) 


so that the expansion point t n+il := t n + y h must be moved into the complex plane. Upon 
substituting y = \ ± i into (3.9), one obtains the complex truncation error 


_ ^ ( t"" k dfn f m\ 

6n+1 ~ 720 5 dy n ' 


(3.12) 


T * 


— h 4 ( f"' - 4 — f") + —h 5 ( f" 

21(3 1 Vn ^ ^ Jn) + ,1 Q o vn 


^/n <•/// 

/ 71 


4 ^ a 


432 " v ' 71 ' % ‘ dy 

The above expression is valid so long as / is analytic in a sufficiently large neighborhood 
of ( t n ,y n ) in the complex plane. 

Equation (3.12) shows that for a real- valued function /, the real truncation error is 
0(h 5 ). This is two orders higher than for any other value of y,. We should note further 
that with y defined as above, (3.6) becomes 
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(3.13a) 


Vn+n 


R e(y n ) + yhf{t n+tl ,y n+fM ) 


(M) 2 

2 


/ (^n+/ji Vn+ii) 


2/n+l 




Vn+ii T A* ^ Un+fi) T 




(3.13b) 


where /i is the complex conjugate of /i. All operations in (3.13) must be performed in 
complex arithmetic, while the solution at each time level is the real part of y n . 

If // is real, one can verify that // = ^ minimizes the truncation error e n+ \ at a value 
of /" + 0(h 4 ), in which case (3.6) becomes 


Jl ft 

y n +i ~ Vn + 2 f(t n +\,yn+\) - -g- f'(t n+ I,y n+ 1 ) (3.14a) 

h h 2 

2/n+l = J/n+j + 2^ <n+ 5 ,2/n+ ^ + -g" /'(*»+*>&.+*)• (3.14b) 

It can also be shown that, corresponding to (3.14), the quadratic expansion (3.1) approx- 
imates y on I n with a maximum error of ^ |/"| + 0(h 4 ). The corresponding error bound 
for both (3.2) and (3.4) is 4^ |/"| + 0(h 4 ). As in the linear case, a centered expansion 
has a maximum error which is one- fourth that of an expansion around either endpoint. 
It should also be observed that, while (3.14) is only a second-order scheme, it provides 
a polynomial approximation to y on I n whose accuracy is one order higher than that of 
( 2 . 12 ). 


Taylor Series of Arbitrary Degree 

Let /( fc )( t , y ) denote the kth total derivative of / with respect to t, where f^(t,y) := 
/(£, y). A Taylor series approximation of degree k to the exact solution y on I n := [t n . £ n +i] 
is 


£ n _j~ 

Vn ~ Vn+ii + ~~ 2/n-b/i) “1“ ? 2/n+M) 


+ ~ f"(t n+fl ,y n+tl ) + ••• + — ^(tn+fi^yn+n), 


(4.1) 


where £ n+#i := £„ + yh and 0 < Re (p) < 1. As shown in the two previous sections, the 
expansion coefficient // defines a family of marching schemes. 

Corresponding to y = 0, one obtains the explicit Taylor series method of order k, 

h? h k 

2/n+l = Vn 4“ hf(t n ,y n ) T ' 2 j "/ {iniVn) + •’• + ^(tn)2/n)- (4-2) 
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Corresponding to /x = 1, one obtains the implicit series method 


Vn+l — Un T Jl f (t n + i . y n + \) f i^n+liVn+l) 

+ ^/"(W UVn+l) + ••• + (-«*-* ( 4 . 3 ) 


When 0 < Re (ji) < 1, one obtains a two-halfstep implicit/explicit marching scheme anal- 
ogous to (2.8) and (3.6), where the implicit and explicit halfsteps are, 

/ ^ \ 2 

Vn+ii ~ Vn + /i /? /(^n-f /x? ?/n+/i) — f 2/n+/x) 


+ J "{ t n + n -. Vn+u) 


+ ... + (4.4a) 


2/n-hl — 2/n+/i H“ (1 AO ^ /(^n+/Ln 2/n+/x) “I" 2 \ ^ i Vn+fi ) 

+ + ••• + 2 izAff/(‘-i)(i„ + „, s „ + „), 

respectively. 

The local error of (4.2) is 

/i fc+1 


,(*) 


/f 1 1 + 


/t fc +2 


h k+s 


J! f(k+ 1 ) , Jl f(k+ 2 ) , 

n+1 (/fe + 1)!'" T (A: + 2 )! Jn (k + 3)\ Jn 


(4.4b) 


(4.5) 


The truncation error for (4.3) and (4.4) can also be derived. Corresponding to (4.3), one 
obtains 


42. = <-u* 1 t£^‘> + I,M 


(* + !)! 


(k + 2)! 


( k + 1) + (k + 2) ^ /(*> 


+ 


h k+3 

(A: + 3)! 


(*+!)(* + 2) 


2 /< l+2 > + (t + l)(i + 3)^/<‘ +1 > 


(A- + 2) (A: + 3) d f' n 

2 /n ' 


+ 0{h k+ 4 ), 


(4.6) 


for all A- > 2. A proof is given in the Appendix. (The error for k — 1 is given by (2.6), and 
only differs from (4.6) in the last term.) 
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The discretization error for (4.4) can be derived in a two-step process. Recognizing 
that (4.4a) is equivalent to (4.3) with a stepsize of fih, one uses (4.6), with h replaced by 
yh, to obtain an expansion for y n+tl . For A odd, one gets 


- » + «/,/„ + M!/: + ... + 


2 ! 


k\ 


(4.7) 


+ ^ h '> k+ 2 ' 


(k + 1)! ' (k + 2)! L ( * + 2) + {k +2) ^f f * k) 


(nh) 


fc+3 


rA 2 + 3A + 4 


df n 


(k + 2) (A + 3) df n 


+ TFT 3j! i 2 f " +2) + (* + l)(» + 3 )^/S‘ +I) + 


+ 0(/i fc+4 ), 


and for A: even, one obtains 


(fc) 

i/n+/i 


+ M/n + 


2 ! 


fn + 


+ 


(/ i ^) fc f(Jfc-l) 

A:! /n 


(/i A) fc+2 

(A + 2)! 




(4.8) 


(y A ) fc + 3 

(A + 3)! 


[ fc (fc^ 3) j(k+2) + (i+lJ^ + S)^/^ 1 ) 


+ 


(A + 2) (A + 3) <9/^ 

2 /n 


+ 0(h k+4 ). 


One then expands f(t n+li , y n + M ) and its total derivatives in (4.4b) in a series about 
(An; Vn) • Powers of (y n+/1 — y n ) that appear are evaluated by way of (4.7) for A odd and 
by (4.8) for A even. It can be shown that y^h agrees with the Taylor series expansion for 
y(t n +i ) through terms of order A. The 0(h k+1 ) term is a function of y, and for A is odd, 
is given by 


ftt+1 h*) 
(A + 1)! Jn 


2 y 


k+l 


+ 


(A + 1) y k (1 - y) 


(4.9) 




(A + 1) A (A — 1) £—2 /i \3 

+ 1 V 


+ ’ ’ • + (A + 1) y (1 — y)^ 
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and for k even, is given by 


hk+i 


(*+!)! 


f {k) [(* + l)/(l-/i) 


(4.10) 


+ 


(k + 1) A: k _ l . 2 , (* + 1) k (k - 1) fc _ 2 ,, ,.\3 , , , n,./i 


2 ! 


A**" 1 (1-/1)- + 


3! 


/i " (1 — /x) + • • • + (k. + 1) /x (1 — /x) 


The sum in brackets in (4.9) and (4.10) can be simplified by making use of the following 
identity, which follows from the binomial theorem. 


/ + */.*->( l-/.) + ^ r ^/- 2 (i-/') 2 


(4.11) 


+ k (k- l)(k - 2) ^_3 (1 _ ^3 + ... + + (!_„)* = 1. 

Equation (4.11) holds for all k > 0. 

Now (4.9) can be written 


h k ' hl 


(Jfe + 1)! 


f {k) < / +1 + k +1 + (Jfc + l)Ai fe (l-/x) 


(4.12) 


+ H k ~' (1 - + ( * + 1 ) 3 , (fe - 1) /- 2 ( 1-^ 3 + + (k + 1 ) H (1 - /«)* 


+ (1-M)‘ +1 -(l-/i) 


1 


It follows from (4.11) that (4.12) equals 

h k +' 

(* + l)! Jn 

Similarly, one rewrites (4.10) as 


f w [/+1 - (1 - /.)*+* + 1 


(4.13) 


h k+l 


(k + 1)! 


f(k) J _^+l + Lfc+l + (A; + l) / x fc (l- / x) 


(4.14) 


+ 


(fcj-l) k k _ x 
2 ! 


„*-* (1 - + (t + 1) 3 t , ( *~ 1) 


/t fc 2 (1 — /() 3 + • ■ • + (k + 1) [i (1 — /i) 2 


+ (1-M)‘ +1 -(!-/•) 


fc+1 
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which equals 


h k+i 


(k + 1)! 


f (k) ^ _ „»+l _ (1 _„)*+! 


(4.15) 


(k) 

Tlie expansion for y n li is then, for k odd or even, 


+ 


Vn + 1 
h k+1 


Vn + h f n + — f' n + •■• + - j /^ fc 


2! Ad 

/^ fc) [l + (-l) fc+1 / +1 - (l-/i) fc+1 l + 0(h k+2 ). 


(4.16) 


( k+l)l 

It follows that the discretization error of the marching scheme (4.4) is 

= [(i - + (-i)V +1 ] (4^1, 

Equation (4.17) is valid for 0 < Re (//) < 1, and for all A; > 1. 

One may also derive the higher order terms for in a manner similar to the above. 
The next two terms in the expansion are 


f {k) + 0{h k+2 ). (4.17) 


(1 - y) k+1 + (*+l)p[(l-/i) fc+1 + (-1)V +1 ] 


f(k+ 1 ) 


+ (— l) fe (A: + 2) fi k+l 


f(k) I / * fc+2 
dy Jn ( (k + 2)< 


(4.18) 


and 


(!-/.)*+■ [(* + i)M + 1] + (* + 2 H t + 1 ) fl 2 [(1 _ /l)t+ . + ( _1)V +1 ] 


+ (-1)* (* + 3) (k+l)^ +2 ^ /<*’+■) 


+ / jNfe (A: + 3) (k + 2) jfc+1 d/4 I h k+3 




dy 


dy 

/< fc) 


f(k+ 2) 

(4.19) 


(A; + 3)!' 


From (4.17) it follows that the optimal expansion coefficients for a Taylor series method 
of degree k are the roots of the equation 


or 


(l-/i) fc+1 + {-l) k n k+1 = 0, 

(4.20) 

-fi) k+1 - fi k+1 = 0 for A; odd, 

(4.21a) 
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and 


(1 - y) k+1 + f.i k+1 =0 for A: even. 


(4.21b) 


Since the left hand sides of (4.21a,b) are polynomials of degree k, each equation has at 
most k distinct roots. One can show that the k roots of (4.21a) are 


1 , . 1 

V = 2 ± l 2 tan 


and the k roots of (4.21b) are 


(m — 1) 7r 
k -f- 1 


m = 1,2, 


k + 1 


(4.22a) 


1 , • 1 . \( m ~ \)^ 
n = - ± i - tan — — 

h 2 2 k + 1 


m = 1, 2, 


k 

2 ' 


(4.22b) 


Equation (4.21a) thus has one real root, y, = while Eq. (4.21b) has no real roots. 

Using (4.17) - (4.20), the discretization error corresponding to all y defined by (4.22) 

is 


o (fc,m) _ 

"71- f 1 


[(i - /‘) t+1 ff +1) + (-i) 1 (* + 2) / +I /!*> 


+ 


h k +3 

(1+3)! L 


(l-x)* +1 [(i + l)/. + l]/f +2) + (-l) fe (fc + 3)(* + l)/ + 2 ^/jS' t+1 > 


, / qk (& + 3 ) ( fc + 2) fc+ i 5/' (fc) 

+ 10 2 [l dy Jn 

Equation (4.23) can be simplified by utilizing, from (4.20), 

(l -^) fc+1 = (-l) fc +V +1 - 


(4.23) 


(4.24) 


One obtains 


o (k,m) 

-n+l 


(■ -i) fc+ V + + 




h k +3 

(k + 3)! 


C(fc +!)[/, { . t+2) - (* + 3) /< t+1) ] 


, f(k+2) _ + 3) (A: + 2) df' n , (fc ) 

+ /n 2 dy Jn 


+ 0(h k+4 ). 


(4.25) 


Now if /r is written in polar form, the expression for can simplified further. 

By defining 6 according to 


9 = ± 


(m — 1) 7r 

k + 1 


for m = 1, 2, • • • , - 1 for A: odd 


(4.26a) 
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and 


e = ± 

(m — i) 7T k „ . 

— - — - — , for m = 1, 2, • • • , — lor k even 

k- 1-1 2 

(4.26b) 

one can write 

p = ^ sec 9 e l ° . 

(4.27) 

One then obtains 

/x fc+1 = ^ sec^^j 

(4.28) 

It follows that for k odd, 

/I \ fc~h 1 

(- sec fl) , 

(4.29a) 

and for k even, 

/I \fc+i 

/.i k+1 = i (— l) m+1 sgn(0) y- secflj 

(4.29b) 


It is significant that f.t k+1 is pure real when k is odd and pure imaginary when k is 
even. An inspection of (4.25) shows that if / is a real-valued function and k is odd, the 
0(h k+2 ) error term is real. On the other hand, if / is real and k is even, the 0(h k+2 ) error 
term is pure imaginary. Consequently, the real truncation error is 0(h k+2 ) when k is odd, 
and 0(h k+3 ) when k is even. One thus obtains to leading order, for all odd k > 1 and for 
/ real, 


a (k,m) 

J n+1 


(-l) t+m (5 


1 x fc+1 },k + 2 

2H (FT2)' L 


4 ‘ +1) - (* + 


(4.30) 


where 6 and m are defined by (4.26a). For even k , one obtains 


Jk,m) 
'n + 1 


(4.31) 


(- 1 ) 


/b+m + 1 


k + 1 


tan 


«l(j 




where 6 and m are defined by (4.26b). 

It is thus clear that, corresponding to the expansion coefficients derived above, a Taylor 
series method can have an accuracy which is one or two orders higher than the classical 
case in which fx — 0. It also turns out that, due to the special form of the truncation error 
in (4.30) and (4.31), the accuracy of the solution can be raised two additional orders by 
local extrapolation. We will return to this point later in the paper. 

One can now simplify the marching scheme (4.4) for all // defined by (4.22). Using 
the fact that 

1 - n = ft, (4.32) 
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(4.4) becomes 


+ 


( j.JL 1 1 ) 2 

Un+n — {Vn) + l 1 h /(^n+M? J/n-fj*) — / (^n+*n 2/n+/z) 

^ - ••• + ••■ + (-l)** 1 n^/“‘ 1) («»+^.»»+|.) ( 4 -33a) 


2/n + l = 2/n-f/x “I" M ^ 2/n-f/i) f Z/n+/i) 




k! 


(4.33b) 


For the special case y = \. one has /i = /i, so that (4.33) can be written 


1 h 2 


Vn+ i = y« + ^/(<n+i^n+i) “ ^ /'(<n+* > */n+i ) 


1 /t 3 „ 

2^ 1}T ^ ^ n + i ’ ^ n +5 ) — 


+ 




(4.34a) 


1 ^ 

IJn+l = Vn + ft f(t n +l , 2/ n -fl) + ”77 / (^n+t>2/n+|) 


1 h 5 , m , . 

2? "5T ^ (^n+i^Ti+i) + 


+ 


2 2 3! 
1 


J_ + (-l)^- 1 — 

2 fe V J 2 k 


h k 


-j^f {k 1} {t n+ i,y n+ i). (4.34b) 


Note that all odd-ordered derivatives have been eliminated in (4.34b). 

Scheme (4.34) is equivalent to a class of Turan type, one-stage implicit Runge-Kutta 
methods [20]. However, unlike implicit Runge-Kutta methods, the implicit step (4.34a) 
only requires the solution of a single equation, instead of a system of k equations. This 
advantage becomes significant when k is much bigger than one. 

From (4.30), the local error of the real marching scheme (4.34) is, for all odd k, 


e 


(k) 
n + 1 


1 h k+ 2 

2^+1 (Jfc + 2)! 


fL k+1) - (* + 2) ^ /<*=> 


+ 0{h k+3 ). 


(4.35a) 


For even k, one obtains from (4.17) 



1 h k+l 

¥ (fc + 1)! 


/f> + 0(h k+2 ). 


(4.35b) 
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One also derives, corresponding to // = the following error bounds on the polynomial 
approximate y n (t ) defined by (4.1). 

k odd: 

I y(t) - Vn(t) | < |/^| + 0(h k + 2 ), t n < t < t n+1 (4.36a) 

k even: 

|?/(0 - Vn(t) | < l/n fc) l + 0(/* fc + 2 ), t n < t < t n+ 1. (4.36b) 

Compared to the corresponding error bounds for \x — 0 and fi = 1, (4.36a) is more accurate 
by a factor of 2 fc+1 , and (4.36b) is more accurate by a factor of 2 k . 


Stability for y' = Ay , Re (A) < 0 


Applying the marching scheme (4.4) to 

y' = Ay, y(t 0 ) = y 0 


one gets 


2/n+l Vn 


1 + (l-n)Xh + (l-y) 2 A 2 f + ••• + (1 — y) fc A fc 
1 - yAh + y 2 A 2 ^f - ••• + + (-l) k y, k X k ~^ 


(5.1) 


(5.2) 


so that the discrete solution is 


Vn 


Vo 


1 + (l-y)Ah + (l-y) 2 A 2 ^ + ••• + (l-y,) k \ k £ 


1 — fi A h + y 2 A 2 


2 ! 


+ • ■ • + (— 1 ) k y, k x k 


hi 

k\ 


(5.3a) 


Equation (5.3a) is valid for all // satisfying 0 < Re (y) < 1. 

Since the exact solution to (5.1), y(t) = y 0 tends to zero as t tends to infinity 

when Re (A) < 0, it follows that the stability criterion is 


1 + (l-g)A/t + (l-y) 2 A 2 f + ••• + (l-y) fc A fc ^- 
1 - t iXh + /x 2 A 2 f + (-l) fc y fc A fe ^- 


We consider the three cases y — 0, y = 1, and Re(y) = for which (5.3) simplifies 
as follows. 
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// = 0 


= Vo ( 1 + A h + A 2 — + • ■■ + A fc — ) 


h k \ n 


Vn 


(5.4a) 


h 2 h k 

1 + A h + A 2 — + • • • + A k — 
2! k\ 


< 1 


(5.4b) 


/' = 1 : 


2/n 


y 0 


.1 - Aft + A*$ - ■■■ + ••• + (-1)»A»£]’ 


1 — X h + A 2 — 


+ ••• + 


> 1 


(5.5a) 


(5.5b) 


Re (p) = \ : 

1 + p Ah + p 2 A 2 ^ • + p k A k ^ 

Vn = y ° 1 - ft A /) + // 2 A 2 ^ -■•• + •■• + (-1)VA A: ^ 

1 + fi A h + /i 2 A 2 • + p k A k < 

1 - pA/i + // 2 A 2 ^f - ••• + ••• + (-l)V fc A fc |^- 


(5.6a) 


(5.6b) 


Degree-1 Taylor Series 

Letting A h = w = u + iv, where u < 0, the stability criterion for Euler’s method is 

|1 + w \ 2 = (1 + M ) 2 + v 2 < 1. (5.7a) 

The above inequality describes the set of points inside a circle of radius one centered at 
the point (-1,0) in the complex Ah plane. For real A, (5.7a) simplifies to 

h < (5.7b) 

— A 

From (5.5b), the stability requirement of the backward Euler method is 

|1 - w\ 2 = (1 - u) 2 + V 2 > 1, (5.8) 

which is clearly satisfied for all negative u. The backward Euler solution thus tends to zero 
and remains stable for any fixed h > 0. Schemes which possess this property for problem 
(5.1) are said to be A-stable [23], [24]. 
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Prom (5.6b), the stability requirement of the centered Euler method is 


1 + [iw 
1 — w 


< 1 , 


where y, = n = i. Since the centered Euler solution 


Vn 


( 1 + Aj V 

= Mrrri) 


(5.9) 


(5.10) 


is identical to that of the trapezoidal rule (which is known to be A-stable), it follows that 
the centered Euler method is also A-stable. (It should be noted that Dahlquist [25] has 
shown that the trapezoidal rule is the most accurate, A-stable scheme from the general 
class of linear multistep methods. The centered Euler method is, of course, a Taylor series 
method and therefore not in the same class. However, both methods give the same solution 
to y' = A y.) 


Degree-2 Taylor Series 

The explicit stability requirement (5.4b) is 


1 + w + 


w 


< 1, 


(5.11a) 


which reduces to 

2u + 2 u 2 + u(u 2 + v 2 ) + l -(u 2 + v 2 ) 2 < 0. (5.11b) 

The above inequality is satisfied for all (a, v) such that 

— 2 < u < 0 (5.12a) 

|t>| < \J 2 \/—u (u + 2) — u{u + 2). (5.12b) 

See Figure 1. 

The fully implicit stability criterion is 


w 


or 


> i, 


i 


(5.13a) 


2 u — 2 u + u(u + v ) — ~^i u + v ) < 0. 


(5.13b) 

Since (5.13b) is satisfied for all u < 0, an implicit Taylor series of degree two is A-stable. 
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Stability requirement (5.6b) is 


1 + (J,w + n 2 ^ 


1 — H W + / L 2 -J- 


< 1. 


(5.14a) 


(5.14b) 


Letting Hi := Im (/x), it can be shown that (5.14a) holds if 

u 2 + 4fiiV + (— + /x 2 ) (tx 2 + u 2 ) < 0. 

Inequality (5.14b) is clearly satisfied for all negative u when /q = 0. For m = ± 
corresponding to (3.12), (5.14b) can be written 


u 


U / V ^ 1 1 \ - 

1 + T + ( T” ± 4 


Vs 


< 0, 


(5.15) 


so that a degree-two series is A-stable for both h = \ and /x = \ ± i -qp. 


Degree- 3 Taylor Series 

Proceeding as above, one derives the following stability constraints for Taylor series 
of degree three. 

/x = 0 : 


H = 1 : 


2u + 2u 2 + + i(u 4 - v 4 ) + j(u 2 + v 2 ) 2 

O O w 

+ ~U (u 2 + v 2 ) 2 + ~ (u 2 + V 2 ) 3 < 0 

O OU 


2u - 2u 2 + ^u 3 - \{u A - v 4 ) - \ {u 2 + v 2 ) 2 

u O vt 

+ i U (u 2 + v 2 ) 2 - (u 2 + V 2 ) 3 < 0 


(5.16) 


(5.17) 


H = \ + i m : 

[2 + 4„.„ + £ + 4rfv 2 + + 5<T + * V ) <T + 


2 

U 


+ g Mi v 


(« 2 + » 2 ) (j + ^ + - 4 ) + "?>] < 


(5.18) 


Because of the complexity of the above expressions, direct analysis is difficult. How- 
ever, one can determine numerically that, corresponding to (5.16), an explicit series of 
degree three has the stability region shown in Figure 2. One can also show that near 
v = ± \/3, (5.17) is not satisfied as it -» 0~, so that a degree-three, fully implicit series is 
not A-stable. See Figure 3. 

Corresponding to (5.18), one considers stability for Hi ~ 0 and Hi = ±\. (See (4.22a).) 
For Hi — 0, one can see directly that (5.18) is always satisfied when u is negative, so that a 
degree-three series is A-stable for H = However, for Hi = one can show that (5.18) 
is not satisfied for all negative u. When /q = + one obtains the instability region 
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y 4 y/— ( v + 2) - (v + 2) 2 < u < 0 


(5.19a) 


-4.5198421 < v < -2, (5.19b) 

and when /x* = — one obtains 

- yWw - 2 - (v - 2) 2 < xx < 0 (5.20a) 

2 < v < 4.5198421. (5.20b) 

See Figure 4. 

Degree-4 Taylor Series 

For Taylor series of degree four, one must rely primarily on numerical means to de- 
termine stability. We obtained the results shown in Figures 5 and 6 for explicit (/x = 0) 
and fully implicit (/i = 1) series, respectively. On the other hand, it is possible to show 
rigorously that, for /x = 4, the degree- four series is A-stable. However, corresponding to 
the two complex pairs /x = ^ ± i ^ tan and /x = | ± i | tan yjj defined by (4.22b), 
the degree-four series is not A-stable. See Figure 7. 

Summary of A-Stable Taylor Series 

To the best of our knowledge, there are no A-stable Taylor series of the form (4.4) 
with degree higher than four. (There are, however, A-stable schemes which use derivatives 
of arbitrarily high order [26], [27].) Table I summarizes the A-stable Taylor series identified 
in this section. 


Taylor Series 

Expansion Coefficient 

Order of Accuracy 

Degree- 1 

/x = 1 (Backward Euler) 

1st Order 


(Centered Euler) 

2nd Order 

Degree-2 

/x = 1 

2nd Order 


/*= h 

2nd Order 


u — i ± i & 
A* 2 31 1 6 

4th Order 

Degree-3 


4th Order 

Degree-4 


4th Order 


Table I A-Stable Taylor series with their orders of accuracy. 
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Numerical Results - Uniform Stepsize 


Our main objectives in this section are to verify the order-of-accuracy results presented 
earlier in the paper and to introduce extrapolation schemes for Taylor series of degree three 
and higher. We also discuss some details concerning the numerical implementation of the 
marching scheme (4.33). For a set of model problems, we consider Problems A1 - A5 from 
Hull et. al. [19], which are shown below. 


Problem 


ODE 

Initial Condition 

Exact Solution 

1. 

y' 

= -y 

2/(0) 

= 1 

y = e ~ t 

2. 

y' 

_ y 3 

2/(0) 

= 1 

1 

2 

v VtT 1 

3. 

y' 

= y COS t 

2/(0) 

= 1 

y = e ,in ‘ 

A 

y' 

= V -{1 -i) 

4 v 20 ' 

2/(0) 

= 1 

20 

4. 

y - l + 19e-*/ 4 

5. 

y' 

y - 1 

2/(0) 

= 4 

r = 4e n / 2 e~ e 


y + 1 


t = r cos 6 
y = r sin 6 






Degree- 1 Taylor Series 

Corresponding to ft = 1 and /i = |, (4.33) becomes the centered Euler method 
(2.12a, b) 

2/n+ i = Vn + j/fWHn+i) 


rt */ x 

Vn+l = Vn+\ "b 2 •U^i+ 2 1 2/n+j)’ 


For nonlinear /, the implicit halfstep (2.12a) must be solved by iteration. One can iterate 
on (2.12a) directly using corrector iteration, 


y 


(0 

n+5 


, h tu (/- 1 )\ 

Vn + 77 Wn+l’Vn+i ) 


( 6 . 1 ) 
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where y 


( 0 ) 

n+5 


= y n . One can alternatively use Newton iteration, 


y 


(0 

n + \ 


y 


(l-l) 

n+i 



( 6 . 2 ) 


For each method, the computational work per step depends on the number of iterations 
M required to solve (2.12a). See Tables II and III. Corresponding results for the implicit 
midpoint rule and trapezoidal rule are also shown for comparison. 

We solved Problems 1 - 5 on the interval [0,20], with stepsize h = 0.5, 0.25, 0.125, .... 
Figure 8 shows the reduction of numerical error with stepsize. The corresponding orders of 
accuracy are shown in Table IV. All calculations were performed in double precision FOR- 
TRAN 77 on a Dell Dimension XPS R400 computer running under the LINUX operating 
system. 


Operation 

Centered Euler 

Implicit Midpoint 

Trapezoidal 

Function 

Evaluations 

M + 1 

M 

M + 1 

Additions and 
Subtractions 

M + 1 

2 M 

2 M 

Multiplications 

M + 1 

2 M 

M 

Divisions 

0 

0 

0 

Table II Computational Work Per Step - Corrector Method 

Centered Euler Implicit Midpoint Trapezoidal 

Operation 

Function 

Evaluations 

2M+1 

2 M 

2 M + 1 

Additions and 
Subtractions 

4 M + 1 

5 M 

5 M 

Multiplications 

2 M + 1 

3 M 

2 M 

Divisions 

M 

M 

M 


Table III Computational Work Per Step - Newton's Method 
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1 

2 

Problem: 

3 

4 

5 

Method 
Centered Euler 

2.00 

1.98 

2.02 

1.99 

2.06 

Trapezoidal Rule 

2.00 

2.01 

2.00 

2.00 

2.02 

Implicit Midpoint 

2.00 

2.00 

1.95 

2.00 

1.98 


Table IV Numerical order of accuracy for Problems 1-5. 


Degree- 2 Taylor Series 

Setting k = 2 and /i = \ ± i \ tan f = \ ± * & in (4.33), one recovers the marching 
scheme (3.13a,b) 


Vn+n — R® (?/n) "h [X h ; Vn+fi) 


Vn + 1 — Vn+iJ 4“ /Z h f (.tfl+ni Vn+n) T 


(»h) 2 „ 


2 

(fib) 2 w 


/ Vn+n) 


f {tn+ni Vn+fi)i 


where the expansion point t n+fJ , t n + fxh is in the complex plane. Taking the “plus” 
sign for fx, the integration path takes the form shown in Figure 9. 

Analogous to (6.1), corrector iteration for (3.13a) is given by 


vS+p = R( - (Vn) + lihf(t n+ p,yZ + * ) ) 


(<-!)> 


(fi hf 


f'(tn+p,yn+, I } ) 


(6.3) 


where y^+n = Re(r/ n )- Similarly, Newton iteration for (3.13a) is 


(Z-l) 

Vn+v 



Letting M denote the number of iterations required to solve (3.13a), the work per 
step is shown in Table V, where all operations must be performed in complex arithmetic. 
Note that the evaluation of in the denominator of (6.4) is not counted as a function 
evaluation, since it is part of the evaluation of f in the numerator. 

Application of the degree-2 Taylor series to Problems 1 - 5 on the interval [0,20] 
produced the error reduction results shown in Figure 10. Table VI shows the corresponding 
orders of accuracy. 
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Corrector Iteration 


Newton Iteration 


Operation 


Function 

2 M + 2 

3 M + 2 

Evaluations 



Additions and 

2 M + 2 

6 M + 2 

Subtractions 



Multiplications 

2 M + 2 

4 M + 2 

Divisions 

0 

M 


Table V Computational work per step for a degree-2 Taylor series with // = | ± i ^ . 


1 


Problem: 

2 3 4 


5 


Expansion 

Coefficient 

+ 4.04 3.96 4.12 3.97 4.23 

Table VI Numerical order of accuracy for Problems 1-5 using a degree-2 Taylor series. 
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Degree-3 Taylor Series 

From (4.22a) and (4.30), the marching scheme (4.33) is fourth-order-accurate when 
k = 3 and 

h = 2 ( 6 ‘ 5 ) 

OI 11 

" = 2 ± *' 2’ (6 ' 0) 

Since // = ^ is real and has the smallest truncation error, it is the preferred expansion 
coefficient. 

The degree-3 Taylor series can be implemented using corrector or Newton iteration, 
analogous to (6.3) and (6.4), respectively. The work per step is shown in Table VII. 

We now show that one can construct an extrapolation scheme using (6.5) - ( 6 . 6 ) which 
is sixth-order-accurate. Let denote the numerical solution at t n+ \ corresponding to 

H — i, and let 

4 3 +i ■= »( ( »+i) - (6-7) 

denote its truncation error. Similarly, let denote the numerical solution at t n + 1 

corresponding top = ^ + i | , and let 


( 3 . 2 ) 
c n + 1 


y{tn+ 1 ) Vn+1 ■ 


( 6 . 8 ) 


(3 i\ 

Then one can see from (4.25) and (4.29a) that for a real-valued function /, er n +i and the 
real part of (^+i are related according to 


Re 



— 4 e 


(3,1) 

n+1 


+ 0(h 7 ). 


(0,9) 


Taking the real part of ( 6 . 8 ) and using (6.9), one gets 

-4e<?;V + 0(h r ) = »(i n+1 ) - Re^i+i’)- (610) 


Multiplying (6.7) by | and (6.10) by | and then adding the resulting equations and 
rearranging, one obtains 

f 4+V + l R*(4+?) = »(*»+i) + o(h 7 ). (o.ii) 

The extrapolated solution 

*- 1 4+i> + 1 r.(*&?) («•«) 

is thus sixth-order-accurate. 

In Figure 11 we show the stability region of the above extrapolation for problem (5.1). 
Figure 12 presents error reduction results for Problems 1-5 for p = 5 , p = \ + 
and extrapolation. The numerical orders of accuracy are shown in Table VIII. 
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Operation 

Corrector Iteration 

Newton Iteration 

Function 

Evaluations 

3 M + 3 

4 M + 3 

Additions and 
Subtractions 

3 M + 3 

8 M + 3 

Multiplications 

3 M + 3 

6 M + 3 

Divisions 

0 

M 

Table VII Computational work per step for 
= 2 ± *5- 

a degree-3 Taylor series with p 



i 

2 

Problem: 

3 

4 

5 

Expansion 

Coefficient 

/* = \ 

3.99 

3.98 

4.01 

3.99 

4.04 

v = i + * § 

4.01 

3.91 

4.02 

3.97 

3.94 

extrapolation 

6.11 

5.89 

6.03 

6.02 

6.20 


or 


Table VIII Numerical order of accuracy for Problems 1-5 using a degree-3 Taylor series. 
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Degree-4 Taylor Series 

From (4.22b), one obtains the expansion coefficients 


and 


h = 




1 .1 7T 

- ± i - tan — 

2 2 10 

(6.13) 

1 1 3 7T 

— i 2 — tan — — , 

2 2 10 

(6.14) 


for which (4.33) is sixth-order-accurate. One can verify from (4.31) that (6.13) has a 
truncation error which is about 47 times smaller than that of (6.14). 

Using a corrector or Newton iteration implementation, the degree-4 Taylor series re- 
quires the work per step shown in Table IX. 

As in the previous section, extrapolation can be used to raise the accuracy of the 
degree-4 Taylor series by two additional orders. Rather than go through the details, we 
simply state the results. 

Let y^+i denote the discrete solution at t n+ 1 when p = | + i \ tan and let y^^\ 
denote the discrete solution when // = \ + i \ tan Define p and q by 


P ■■= tan ^ 



and 



Then the extrapolated solution 


satisfies 


(4,1) , (4,2) 

(4,e) On+ 1 + PVn± 1 

Vn+1 ‘ p + q 

y { n£ = y(t n+ 1 ) + 0(h Q ) 


(6.15a) 


(6.15b) 


(6.16) 

(6.17) 


so that y^+i i s eighth-order-accurate. 

The stability region of the above extrapolation for problem (5.1) is shown in Figure 
13. Figure 14 shows reduction of error with stepsize for Problems 1-5 using p = 
\ + i \ tan p = \ + i \ tan and extrapolation. Table X shows the corresponding 
orders of accuracy. 
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Operation 

Corrector Iteration 

Newton Iteration 

Function 

Evaluations 

4 M + 4 

5 M + 4 

Additions and 
Subtractions 

4 M -|- 4 

10 M + 4 

Multiplications 

4 M + 4 

8 M + 4 

Divisions 

0 

M 


Table IX Computational work per step for a degree-4 Taylor series with 
/x = \ ± i \ tan ^ or /x = ^ ± i ^ tan . 


Problem: 



1 

2 

3 

4 

5 

Expansion 

Coefficient 






/x = | + * J tan ^ 

6.00 

5.94 

6.01 

5.99 

5.89 

/X = j + *2 ^ an To” 

6.00 

5.80 

5.92 

6.00 

6.27 

extrapolation 

7.98 

7.62 

8.01 

7.97 

7.89 


Table X Numerical order of accuracy for Problems 1-5 using a degree-4 Taylor series. 
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Extrapolation for Taylor Series of Arbitrary Degree 


According to our work in the two previous sections, a degree-3 Taylor series with ex- 
trapolation provides a sixth-order method, and a degree-4 Taylor series with extrapolation 
provides an eighth-order method. In general, if k is odd, a degree-fc Taylor series with 
extrapolation provides a (k + 3)rd-order method, and if k is even, a degree-fc Taylor series 
with extrapolation provides a (k + 4)th-order method. The purpose of this section is to 
generalize these results. 

Let denote the extrapolated solution of degree k at t n+i . Let yl^ and 

denote the regular solution of degree k at t n+ 1 corresponding to p = p\ and p = pz, 
respectively, where p\ and /Z 2 are defined as follows. For k odd, let 





(6.18a) 


and 


For k even let 


and 


11 7T 

«2 = - + i - tan 

h 2 2 k + 1 


1,1 - 5 + ‘ \ tan JotTT) 

1 1 37T 

^2 = 2+ a 2 tan 2(*TT)- 


The extrapolated solution for all k is then 


(6.18b) 

(6.19a) 

(6.19b) 


(k,e) 

Vn+l 


(k, 1) , (fc,2) 

QVn+1 + PVn+l 


p + q 

where p and q are defined below. For all odd k > 3, 


( 6 . 20 ) 


p = 1 


(6.21a) 


and 


For all even k > 4, 


and 


7T \ fc+1 


= ( SeC iTl) 


p = tan 


7T 


7 r 


sec 


2(k+l) L 2(Jfc + l)J 


fc+i 


3 7r r 3 7r ifc-i- 1 
tan rv I sec 


2{k+l) r^2(fc + l)] 


(6.21b) 

(6.22a) 

(6.22b) 


Using the above formulas, we implemented extrapolation schemes of degree five and 
six. Figure 15 shows the corresponding stability regions for problem (5.1). Application 
to Problem 1 produced the results shown in Figures 16. a, b, where the numerical orders of 
accuracy are 8.01 and 9.87, respectively. 
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Numerical Results - Variable Stepsize 


The purpose of this section is to briefly introduce a variable step approach that can be 
implemented using an implicit Taylor series method. Here we present results for Problems 
1-5, and in the next two sections we consider singular equations and stiff systems. 

To develop a variable step scheme, one must have a method for estimating the local 
error and a procedure for controlling the stepsize. For an implicit series of degree three 
or higher, the local error can be obtained by extrapolation. For a degree-two series, the 
error can be obtained by integrating twice at each step - once with (3.13) and once with 
(3.14). The difference between the two solutions is an 0(h 3 ) error estimate. Similarly, one 
integrates with (2.12) and (2.5) when using a series of first degree. 

The stepsize may be adjusted as follows. Let e n denote the error at step n, and let r 
be the error tolerance. Then if e n is 0(h p ) and h n is a successful stepsize, an estimate for 
h n +i is 

h n +i = thn (-Y (7.1) 

\ e n / 

where £ is an adjustable parameter less than one [28]. (See [29] for a discussion of stepsize 
control for long Taylor series methods.) 

In what follows, we present numerical results from four variable step schemes, which 
we identify as follows. Scheme 1 is the degree-1 Taylor series (2.12) combined with (2.5). 
Scheme 2 is the degree-2 series (3.13) combined with (3.14). Schemes 3 and 4 are the 
degree-3 and -4 series with extrapolation, respectively. See Table XI. 


Scheme 

1 

2 

3 

4 


Taylor Series Order of Accuracy q_ 

Degree 1 2nd 2 

Degree 2 4th 3 

Degree 3 6th 5 

Degree 4 8th 7 


Table XI Variable step schemes with their orders of accuracy and p values. 


We integrated Problems 1 - 5 on the interval [0,20]. For Schemes 2-4, the initial 
stepsize was determined from < t, where k is the degree of the corresponding 

Taylor series and the total derivative f^ k ~^ was evaluated near t 0 but not at t 0 . For 
Scheme 1, we used \ A 2 < r, where A = [30]. Subsequent stepsizes were determined 

from (7.1), subject to < hf ac t, where hf ac t = 2.0. We took £ to be 0.9 for Schemes 2 - 
4 and 0.7 for Scheme 1. Unsuccessful stepsizes h' n+1 were reduced using h n +i := vh' n+l , 
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where v — 0.7. At the initial step we took v = 0.1. All parameters identified above were 
maintained at their stated values for all calculations. 

For each problem, we varied the tolerance from 10 -2 to 10 -12 , and determined the 
maximum absolute error and number of integration steps per tolerance. The results are 
shown in Tables XII - XVI. The average number of repeat steps per tolerance is stated at 
the bottom of each table. 


Scheme: 



1 

2 

3 

4 

Tolerance 





0.10E-01 

0.0740 

0.0437 

0.0291 

0.0150 

0.10E-02 

0.0789 

0.0453 

0.0328 

0.0209 

0.10E-03 

0.0807 

0.0447 

0.0331 

0.0212 

0.10E-04 

0.0814 

0.0454 

0.0327 

0.0212 

0.10E-05 

0.0816 

0.0454 

0.0329 

0.0210 

0.10E-06 

0.0816 

0.0452 

0.0327 

0.0213 

0.10E-07 

0.0817 

0.0450 

0.0329 

0.0208 

0.10E-08 

0.0817 

0.0439 

0.0331 

0.0200 

0.10E-09 

0.0817 

0.0355 

0.0297 

0.0196 

0.10E-10 

0.0817 

0.0229 

0.0231 

0.0171 

0.10E-11 

0.1690 

0.0128 

0.0166 

0.0147 

Table XILa 

Maximum error in units 

i of the error tolerance for Problem 1 



Scheme: 




1 

2 

3 

4 

Tolerance 





0.10E-01 

24 

11 

8 

7 

0.10E-02 

68 

18 

11 

9 

0.10E-03 

205 

32 

15 

11 

0.10E-04 

642 

61 

21 

14 

0.10E-05 

2023 

124 

31 

17 

0.10E-06 

6391 

258 

46 

23 

0.10E-07 

20204 

545 

70 

30 

0.10E-08 

63886 

1164 

108 

39 

0.10E-09 

202023 

2497 

167 

53 

0.10E-10 

638849 

5368 

260 

71 

0.10E-11 

2020214 

11553 

407 

96 


Table XILb Number of integration steps per tolerance for Problem 1. Average number 
of repeat steps is 0.00 for all four schemes. 
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Scheme: 



1 

2 

3 

4 

Tolerance 





0.10E-01 

0.0447 

0.0179 

0.0129 

0.0026 

0.10E-02 

0.0509 

0.0195 

0.0330 

0.0127 

0.10E-03 

0.0530 

0.0149 

0.0447 

0.0220 

0.10E-04 

0.0536 

0.0090 

0.0491 

0.0312 

0.10E-05 

0.0538 

0.0047 

0.0463 

0.0389 

0.10E-06 

0.0539 

0.0023 

0.0391 

0.0439 

0.10E-07 

0.0539 

0.0011 

0.0300 

0.0453 

0.10E-08 

0.0539 

0.0005 

0.0217 

0.0430 

0.10E-09 

0.0539 

0.0003 

0.0149 

0.0384 

0.10E-10 

0.0777 

0.0004 

0.0100 

0.0325 

0.10E-11 

0.2639 

0.0043 

0.0067 

0.0267 


Table XHI.a Maximum error in units of the error tolerance for Problem 2. 


Scheme: 


1 


Tolerance 


0.10E-01 

22 

0.10E-02 

62 

0.10E-03 

189 

0.10E-04 

593 

0.10E-05 

1868 

0.10E-06 

5899 

0.10E-07 

18649 

0.10E-08 

58968 

0.10E-09 

186466 

0.10E-10 

589650 

0.10E-11 

1864632 


2 

3 

4 

10 

8 

7 

17 

10 

8 

31 

13 

10 

60 

18 

13 

121 

25 

16 

252 

37 

20 

535 

54 

26 

1143 

82 

34 

2453 

125 

45 

5275 

194 

60 

11352 

303 

81 


Table XIII. b Number of integration steps per tolerance for Problem 2. Average number 
of repeat steps is 0.00 for all four schemes. 
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Scheme: 


1 


Tolerance 


0.10E-01 

0.3811 

0.10E-02 

0.4194 

0.10E-03 

0.5615 

0.10E-04 

0.6472 

0.10E-05 

0.9096 

0.10E-06 

0.8700 

0.10E-07 

0.9516 

0.10E-08 

1.0210 

0.10E-09 

1.1632 

0.10E-10 

1.5815 

0.10E-11 

6.7590 


2 

3 

4 

0.1431 

2.7849 

1.4753 

0.1012 

1.6676 

6.0293 

0.0338 

0.4951 

2.9164 

0.0394 

0.3809 

0.9818 

0.0187 

0.0843 

0.6711 

0.0112 

0.1057 

2.0509 

0.0072 

0.1928 

0.2164 

0.0035 

0.0881 

0.3157 

0.0019 

0.0953 

0.3654 

0.0139 

0.0205 

0.2030 

0.2318 

0.0355 

0.1470 


Table XIV. a Maximum error in units of the error tolerance for Problem 3 


Scheme: 


1 


Tolerance 


0.10E-01 

189 

0.10E-02 

583 

0.10E-03 

1833 

0.10E-04 

5785 

0.10E-05 

18282 

0.10E-06 

57802 

0.10E-07 

182774 

0.10E-08 

577969 

0.10E-09 

1827690 

0.10E-10 

5779651 

0.10E-11 

18276843 


2 

3 

4 

48 

24 

18 

91 

35 

24 

182 

54 

32 

374 

81 

45 

786 

124 

59 

1676 

189 

78 

3587 

293 

109 

7706 

457 

148 

16580 

718 

200 

35699 

1129 

275 

76885 

1784 

376 


Table XIV. b Number of integration steps per tolerance for Problem 3. Average number 
of repeat steps is 13.27, 39.55, 21.55, and 18.64, respectively. 
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Scheme: 



1 

2 

3 

4 

Tolerance 





0.10E-01 

0.5808 

0.0845 

0.1886 

0.1124 

0.10E-02 

0.6372 

0.0447 

0.0563 

0.2916 

0.10E-03 

0.7465 

0.0358 

0.0960 

0.1249 

0.10E-04 

0.7706 

0.0253 

0.1353 

0.1802 

0.10E-05 

0.8511 

0.0176 

0.1521 

0.1105 

0.10E-06 

0.8976 

0.0132 

0.1325 

0.0940 

0.10E-07 

0.9536 

0.0070 

0.1009 

0.0655 

0.10E-08 

1.0350 

0.0049 

0.0670 

0.0680 

0.10E-09 

1.0712 

0.0031 

0.0706 

0.0458 

0.10E-10 

1.3559 

0.0064 

0.0375 

0.0276 

0.10E-11 

6.3949 

0.1315 

0.0426 

0.0249 


Table XV. a Maximum error in units of the error tolerance for Problem 4. 


Scheme: 


1 2 3 4 

Tolerance 


0.10E-01 

61 

15 

9 

8 

0.10E-02 

188 

26 

12 

9 

0.10E-03 

587 

48 

16 

11 

0.10E-04 

1851 

94 

22 

13 

0.10E-05 

5845 

193 

31 

17 

0.10E-06 

18476 

406 

45 

21 

0.10E-07 

58421 

864 

68 

27 

0.10E-08 

184736 

1849 

102 

35 

0.10E-09 

584179 

3973 

157 

46 

0.10E-10 

1847330 

8546 

244 

61 

0.10E-11 

5841766 

18398 

382 

81 


Table XV. b Number of integration steps per tolerance for Problem 4. Average number 
of repeat steps is 1.82, 5.64, 2.00, and 2.45, respectively. 
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Scheme: 



1 

2 

3 

4 

Tolerance 





0.10E-01 

0.8047 

0.4356 

0.2230 

0.0720 

0.10E-02 

0.8404 

0.3574 

0.4007 

0.1866 

0.10E-03 

0.8525 

0.2637 

0.4884 

0.3063 

0.10E-04 

0.8564 

0.1598 

0.5048 

0.4553 

0.10E-05 

0.8576 

0.0965 

0.4550 

0.5675 

0.10E-06 

0.8580 

0.0512 

0.3541 

0.5609 

0.10E-07 

0.8582 

0.0293 

0.2658 

0.5781 

0.10E-08 

0.8584 

0.0154 

0.1848 

0.5165 

0.10E-09 

0.9120 

0.0109 

0.1254 

0.4395 

0.10E-10 

0.6662 

0.0266 

0.0828 

0.3670 

0.10E-11 

12.1506 

0.1926 

0.0400 

0.3304 


Table XVI. a Maximum error in units of the error tolerance for Problem 5. 


Scheme: 


1 2 3 4 

Tolerance 


0.10E-01 

66 

13 

8 

7 

0.10E-02 

203 

23 

11 

8 

0.10E-03 

635 

41 

15 

10 

0.10E-04 

2003 

82 

21 

13 

0.10E-05 

6328 

169 

30 

16 

0.10E-06 

20005 

355 

45 

21 

0.10E-07 

63254 

757 

68 

28 

0.10E-08 

200021 

1620 

104 

36 

0.10E-09 

632515 

3479 

161 

49 

0.10E-10 

2000183 

7484 

251 

65 

0.10E-11 

6325129 

16113 

393 

88 


Table XVI.b Number of integration steps per tolerance for Problem 5. Average number 
of repeat steps is 0.00, 3.36, 0.00, and 0.00, respectively. 
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Numerical Results - ODE’s with a Singular Point 


Three ODE's with a regular singularity at the origin are Bessel’s equation, 

x 2 y" + xy' + (x 2 - k 2 ) y = 0 (8.1) 

Coulomb’s equation, 

x 2 y" + [— L (L + 1) — 2 rjx + x 2 ]y = 0 (8.2) 

and the hypergeometric equation (Gauss’s equation) 

x(l—x)y" + [7 — (a + fl + 1) x\ y' — a fly = 0. (8-3) 

The above equations are typical of a large class of linear, second-order equations from 
mathematical physics with a regular singular point. A general approach for integrating 
such equations has been presented by Holubec and Stauffer [31]. Haftel et. al. [32] presented 
a similar approach for coupled systems. We show here that implicit Taylor series methods 
can directly integrate such problems as a special case. We consider the following set of 
model problems. 


Problem 

ODE 

Initial Condition 

Exact Solution 

6. a 

(8.1) with k = 0 

2/(0) 

= 1 

Jo 



2/(0) 

= 0 


6.b 

(8.1) with k = 1 

2/(0) 

= 0 

J 1 



2/(0) 

= 1/2 


7. 

(8.2) with L = 0, 

2/(0) 

= 0 

Coulomb 


rj = 1 

1 

II 

o' 

f \J tx / sinh7r 

Wavefunction 

8. a 

(8.3) with a — 0.4, 

2/(0) 

= 1 

Hypergeometric 


fl = 1.0, 7 = 0.5 

2/(0) 

= 0.8 

Function 

8.b 

(8.3) with a = 0.5, 

2/(0) 

= 1 

Hypergeometric 


fl - 0.5, 7 = 1.5 

2/(0) 

= VO 

Function 


We first looked at order of accuracy near the singular point. Each problem was 
reformulated as a system of two first order equations in the usual way. Problems 6 and 7 
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were integrated on [0,1], and Problem 8 on [0, ^]. Error reduction results are plotted in 
Figures 17 - 19, and the corresponding orders of accuracy are summarized in Table XVII. 

It is evident that a significant loss of accuracy occurred on Problems 6.b and 7. The 
degree-4 Taylor series lost from one to four orders of accuracy, while the degree-2 and -3 
series lost from one to two orders. However, it should be noted that the orders of accuracy 
above are as high or higher than an explicit Taylor series of the same degree would be 
for a nonsingular problem. We also observe from Figures 17 - 19 that without exception 
the extrapolated solution for the degree-3 and -4 Taylor series is more accurate than the 
non-extrapolated solutions. 

To look at variable step performance, we used Schemes 1 - 4 to solve Problems 6.a,b 
on the interval [0,20]. Each scheme was used as described in Section 7, with appropriate 
modifications for a system. The results shown in Tables XVIII - IXX are indeed comparable 
to those in Tables XII - XVI for the nonsingular Problems 1-5. 



6. a 

6.b 

Problem: 

7 

8. a .. 

8.b 

Tavlor Series 






Degree- 1 






v = \ 

2.00 

1.89 

1.85 

2.01 

1.92 

Degree- 2 






u - I + 

3.93 

2.00 

3.02 

3.82 

4.07 

Degree-3 






h = \ 

3.99 

3.89 

3.82 

4.01 

3.88 

p = i + 

3.99 

3.89 

3.81 

3.84 

3.77 

extrapolation 

5.83 

4.47 

4.21 

5.57 

5.77 

Degree-4 






= 2 + * 2 ^ an To 

5.92 

4.01 

4.99 

5.84 

5.75 

/r 2 T ® 2 ^ an TV 

5.92 

4.10 

5.02 

5.70 

5.62 

extrapolation 

5.83 

3.94 

4.83 

7.20 

7.48 


Table XVII Order of accuracy near a regular singular point. 
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Scheme: 



1 

2 

3 

4 

Tolerance 





0.10E-01 

1.0965 

1.7756 

0.4928 

1.1900 

0.10E-02 

1.0755 

0.5731 

0.4334 

1.0079 

0.10E-03 

1.0736 

0.2519 

0.2957 

0.4966 

0.10E-04 

1.0733 

0.1157 

0.2076 

0.3070 

0.10E-05 

1.0732 

0.0536 

0.1335 

0.2350 

0.10E-06 

1.0732 

0.0249 

0.0849 

0.1643 

0.10E-07 

1.0732 

0.0116 

0.0543 

0.1179 

0.10E-08 

1.0733 

0.0054 

0.0340 

0.0878 

0.10E-09 

1.0744 

0.0024 

0.0216 

0.0624 

0.10E-10 

1.1015 

0.0012 

0.0130 

0.0622 

0.10E-11 

1.3933 

0.0212 

0.0106 

0.0892 


Table XVIII. a Maximum error in units of the error tolerance for Problem 6. a. 


Scheme: 


1 2 3 4 

Tolerance 


0.10E-01 

114 

24 

15 

13 

0.10E-02 

332 

51 

23 

18 

0.10E-03 

1021 

107 

34 

22 

0.10E-04 

3198 

226 

52 

28 

0.10E-05 

10083 

481 

81 

37 

0.10E-06 

31854 

1032 

125 

49 

0.10E-07 

100701 

2217 

196 

66 

0.10E-08 

318415 

4769 

308 

89 

0.10E-09 

1006888 

10269 

484 

120 

0.10E-10 

3184029 

22115 

764 

164 

0.10E-11 

10068755 

47639 

1206 

225 


Table XVIII. b Number of integration steps per tolerance for Problem 6. a. Average 
number of repeat steps is 0.00, 0.09, 0.64, and 0.45, respectively. 
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Scheme: 



1 

2 

3 

4 

Tolerance 





0.10E-01 

1.1825 

2.0417 

0.5799 

3.3158 

0.10E-02 

1.1569 

0.6318 

0.4310 

0.8488 

0.10E-03 

1.1535 

0.2763 

0.3510 

0.5983 

0.10E-04 

1.1521 

0.1262 

0.2365 

0.3736 

0.10E-05 

1.1511 

0.1003 

0.1515 

0.2587 

0.10E-06 

1.1505 

0.1007 

0.0958 

0.1845 

0.10E-07 

1.1502 

0.0993 

0.0604 

0.1307 

0.10E-08 

1.1501 

0.0941 

0.0378 

0.0940 

0.10E-09 

1.1492 

0.0832 

0.0237 

0.0736 

0.10E-10 

1.1454 

0.0720 

0.0172 

0.2112 

0.10E-11 

1.3087 

0.0774 

0.0201 

0.5614 


Table IXX.a Maximum error in units of the error tolerance for Problem 6.b. 


Scheme: 


1 2 3 4 

Tolerance 


0.10E-01 

113 

24 

18 

19 

0.10E-02 

326 

50 

25 

23 

0.10E-03 

1000 

106 

37 

27 

0.10E-04 

3131 

224 

55 

34 

0.10E-05 

9869 

477 

83 

42 

0.10E-06 

31177 

1022 

127 

54 

0.10E-07 

98557 

2194 

197 

71 

0.10E-08 

311633 

4718 

308 

93 

0.10E-09 

985438 

10152 

484 

125 

0.10E-10 

3116197 

21855 

761 

169 

0.10E-11 

9854247 

47063 

1200 

229 


Table IXX.b Number of integration steps per tolerance for Problem 6.b. Average number 
of repeat steps is 0.00, 0.09, 0.09, and 0.00, respectively. 
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Numerical Results - Stiff Equations 


In this section we briefly consider stiff equations. For a set of model problems we use 
Problems 9-13 below T . 

Problem 9: (Enright and Pryce [33], Problem A3) 

Vi = -10 4 yi + 100 2/2 - 10 2/3 + 2/4 
2/2 = — 10 3 2/2 + 10 y 3 ~ 10 2/4 
2/3 = -2/3 + 10 y 4 
2/4 = — 0.1 2/4 

yi(0) = l, 2/2(0) = 1, 2/3(0) = 1, 2 / 4 ( 0 ) = 1, t end = 20 


Problem 10: ([33], Problem Bl) 

y[ = -2/1 + 2/2 
2/2 - — 100 2/1 - 2/2 
?/3 = — 100 2/3 + y 4 
Va = - 10 4 2/3 - 100 y 4 

2/l(0) = 1, 2/2(0) =0, 2/3(0) = 1, 2/4(0) =0, £end = 20 

Problem 11: ([33], Problem D4) 

y[ = -0.0131/1 - 1000 2 / 1 2/3 

y' 2 = - 2500 2/2 2/3 

2/3 = - 0.013 yi - 1000 2/1 2/3 - 2500 ?/2 2/3 

2/l(0) = 1, 2/2(0) = l, 2/3(0) =0, tend, — 30 
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Problem 12: (Robertson [34]) 


y[ = — 0.04 ?/! + 10 4 t /2 2/3 

y ' 2 = 0.04 2/1 - 10 4 2/2 2/3 - 3 • 10 7 y\ 

Vz = 3 • 10 2/2 

2/1 (0) = 1, 2 / 2 ( 0 ) = 0, 2 / 3 ( 0 ) = 0, iend = 40 


Problem 13: ([33], Problem F3) 

y [ = -10 7 2/2 2/1 + 10 2/3 

2/2 = — 10 7 2/2 2/1 - 10 7 2/2 2/5 + 10 2/3 + 10 2/4 

2/3 = 10 7 2/2 2/1 - 1.001 -10 4 2/3 + 10 -3 2/4 

2/4 = 10 4 2/3 - 10 . 001 2/4 + 10 7 2/2 2/5 

2/5 = 10 2/4 - 10 7 2/2 2/5 


2/i(0) = 4 ■ 10~ 6 , 2/2(0) = 1 • 10 “ 6 , 


2/3(0) = 0, 2/4(0) = 0, 2/5(0) = 0 


tend = 100 


We integrated each system using Schemes 1-4, along with a trapezoidal extrapolation 
scheme which we now introduce. 

Let 2 /„ denote the trapezoidal solution (2.15), and let y ^ denote the degree-2 Taylor 
series solution (3.14) with /x = 5 . One can show that the extrapolated solution 


Vn+l : = 



( 9 . 7 ) 
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is A-stable and fourth-order-accurate with truncation error 


V{tn+l) ~ Vn+l = ~ 


h 5 

2880 L 


df n 


d 2 fn , d 2 f n 


f/ltl I r w -> n fill , on ( Jn _L 

fn + 5 lJ-/n + 20 V oTHT + 


dij 


dtdy dy 2 


/n) /, 


(9.8) 


A variable stepsize for (9.7) can be obtained using (7.1) together with the 0(h 3 ) error 
estimate e n = |y£ +1 — The above extrapolation with variable step implementation 

will be called Scheme 2T. 

In integrating Problems 9-13, the parameters £, and v, introduced in Section 
7, were assigned default values of 0.9, 2.0, and 0.7, respectively, except for Scheme 1 where 
£ was set to 0.7. Any deviation from these values is noted at the bottom of Tables XX 
- XXIV. The local error estimate e n was taken to be the largest absolute error in any 
solution component. The initial stepsize w T as determined as in Section 7 (with appropriate 
modifications for a system). 

Tables XX - XXI demonstrate some degradation in the ability of Schemes 1 - 4 to 
meet the error tolerance. Whether this holds in general for stiff problems requires further 
investigation. 

One also observes that the eighth-order-accurate Scheme 4 frequently required more 
steps than some of the lower order schemes. This contrasts with the nonstiff results pre- 
sented earlier, and illustrates the well known need for stiff solvers to be able to vary their 
order. 


Scheme: 



1 

2 

2 T 

3 

4 

Tolerance 






0.10E-01 

0.1512 

0.0485 

0.0096 

0.0326 

0.3952 

0.10E-02 

0.2270 

0.0628 

0.0129 

0.0332 

1.0237 

0.10E-03 

0.3221 

0.1356 

0.0164 

0.0339 

2.1679 

0.10E-04 

0.3287 

0.1016 

0.0229 

0.0393 

2.8883 

0.10E-05 

0.3292 

0.0450 

0.0155 

0.0488 

6.2514 

0.10E-06 

0.3293 

0.0453 

0.0100 

0.0695 

4.1981 

0.10E-07 

0.3293 

0.0455 

0.0095 

0.4929 

6.7522 

0.10E-08 

0.3295 

0.0443 

0.0095 

0.5759 

5.6286 

0.10E-09 

0.3265 

0.0327 

0.0086 

1.9309 

6.1883 

0.10E-10 

0.4206 

0.0176 

0.0061 

5.1021 

9.4324 

0.10E-11 

1.4309 

0.1217 

0.3659 

5.7487 

5.8660 


Table XX. a Maximum error in units of the error tolerance for Problem 9. 
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Scheme: 



1 

2 

2 T 

3 

4 

Tolerance 






0.10E-01 

121 

34 

34 

24 

31 

0.10E-02 

369 

60 

60 

32 

34 

0.10E-03 

1157 

116 

116 

44 

41 

0.10E-04 

3654 

236 

236 

65 

53 

0.10E-05 

11551 

496 

496 

96 

72 

0.10E-06 

36525 

1055 

1055 

145 

105 

0.10E-07 

115498 

2259 

2259 

224 

163 

0.10E-08 

365234 

4851 

4851 

349 

249 

0.10E-09 

1154967 

10438 

10438 

547 

406 

0.10E-10 

3652325 

22473 

22473 

883 

662 

0.10E-11 

11549662 

48402 

48402 

1450 

1119 


Table XX. b Number of integration steps per error tolerance for Problem 9. Average 
number of repeat steps is 0.18, 0.91, 1.18, 5.82, and 13.00, respectively. Changes to default 
values: Scheme 3, v = 0.5; Scheme 4, hf act = 1.5, v = 0.5. 





Scheme: 




1 

2 

2 T 

3 

4 

Tolerance 






0.10E-01 

1.5593 

2.8866 

0.4222 

0.7640 

7.1984 

0.10E-02 

2.0085 

3.0205 

0.6115 

0.8367 

3.7537 

0.10E-03 

2.0902 

2.1121 

0.6816 

0.5848 

3.2502 

0.10E-04 

2.1074 

3.2174 

0.4245 

0.9580 

11.8994 

0.10E-05 

2.1121 

3.3088 

0.7441 

0.6915 

4.3284 

0.10E-06 

2.1134 

2.2518 

0.7997 

0.5903 

3.3849 

0.10E-07 

2.1141 

2.8611 

0.4129 

0.8654 

2.0479 

0.10E-08 

2.1142 

0.9528 

0.2677 

0.3111 

2.3132 

0.10E-09 

2.1166 

0.4421 

0.1545 

0.2352 

0.4552 

0.10E-10 

2.1143 

0.2346 

0.1879 

0.2052 

0.2391 

0.10E-11 

17.8417 

0.4707 

2.7374 

0.1712 

0.1882 

Table XXI. a 

Maximum error in 

units of the 

error tolerance 

for Problem 10. 
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Scheme: 



1 

2 

2 T 

3 

4 

Tolerance 






0.10E-01 

712 

135 

142 

73 

44 

0.10E-02 

2219 

299 

304 

118 

70 

0.10E-03 

6992 

642 

648 

192 

105 

0.10E-04 

22082 

1380 

1389 

307 

147 

0.10E-05 

69805 

2969 

2978 

486 

211 

0.10E-06 

220720 

6395 

6397 

766 

297 

0.10E-07 

697952 

13771 

13775 

1208 

413 

0.10E-08 

2207094 

29655 

29656 

1905 

571 

0.10E-09 

6979416 

63864 

63864 

3008 

792 

0.10E-10 

22070843 

137568 

137569 

4757 

1092 

0.10E-11 

69794029 

296360 

296361 

7528 

1512 


Table XXI.b Number of integration steps per error tolerance for Problem 10. Average 
number of repeat steps is 19.45, 26.64, 26.73, 19.18, and 10.00, respectively. Changes to 
default values: Scheme 1, v = 0.5; Scheme 2, £ = 0.8; Scheme 2T, £ = 0.8; Scheme 3, £ 
= 0.8; Scheme 4, £ = 0.8, v — 0.5. 





Scheme: 




1 

2 

2 T 

3 

4 

Tolerance 






0.10E-01 

23 

13 

13 

15 

40 

0.10E-02 

31 

15 

15 

16 

44 

0.10E-03 

56 

17 

17 

17 

46 

0.10E-04 

135 

21 

21 

18 

44 

0.10E-05 

387 

29 

29 

19 

46 

0.10E-06 

1187 

46 

46 

22 

48 

0.10E-07 

3720 

83 

83 

26 

51 

0.10E-08 

11734 

164 

164 

32 

54 

0.10E-09 

37081 

337 

337 

43 

77 

0.10E-10 

117237 

710 

710 

60 

117 

0.10E-11 

370738 

1515 

1515 

89 

189 


Table XXII Number of integration steps per error tolerance for Problem 11. Average 
number of repeat steps is 0.00, 1.45, 1.55, 0.09, and 10.36, respectively. Changes to default 
values: Scheme 3, v = 0.5; Scheme 4, = 0.8, hf act = 1.5, v = 0.5. 
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Scheme: 



1 

2 

2 T 

3 

4 . 

Tolerance 






0.10E-01 

29 

15 

26 

41 

85 

0.10E-02 

48 

91 

22 

43 

79 

0.10E-03 

109 

715 

25 

45 

90 

0.10E-04 

303 

41 

41 

47 

73 

0.10E-05 

920 

75 

77 

47 

72 

0.10E-06 

2874 

153 

154 

49 

77 

0.10E-07 

9054 

319 

321 

60 

74 

0.10E-08 

28599 

676 

679 

78 

78 

0.10E-09 

9040G 

1448 

1450 

112 

84 

0.10E-10 

285857 

3108 

3111 

164 

115 

0.10E-11 

903936 

6686 

6688 

249 

168 

Table XXIII 

Number of integration steps per 

error tolerance for Problem 12. Average 

number of repeat steps is 0.00, 3.91, 

2.82, 5.55, and 28.18, 

respectively. Changes to default 

values: Scheme 2, v = 0.5; Scheme 2T, v — 0.5; 

Scheme 3, £ = 0.8, hf act = 

1.5, v — 0.5; 

Scheme 4, £ = 

0.8, /ifact — 1-0, V — 

0.5. 








Scheme: 




1 

2 

2 T 

3 

4 

Tolerance 






0.10E-01 

25 

11 

11 

14 

32 

0.10E-02 

26 

13 

13 

15 

30 

0.10E-03 

28 

15 

15 

17 

31 

0.10E-04 

30 

16 

16 

18 

36 

0.10E-05 

32 

18 

18 

19 

34 

0.10E-06 

41 

21 

20 

20 

36 

0.10E-07 

77 

28 

27 

23 

37 

0.10E-08 

192 

43 

43 

29 

39 

0.10E-09 

561 

78 

77 

38 

43 

0.10E-10 

1729 

151 

150 

53 

50 

0.10E-11 

5425 

308 

307 

78 

61 


Table XXIV Number of integration steps per error tolerance for Problem 13. Average 
number of repeat steps is 0.00, 0.00, 0.00, 0.00, and 4.09, respectively. Changes to default 
values: Scheme 4, £ = 0.8, hf act = 1.5, u = 0.5. 
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Discussion and Conclusion 


In this paper we have presented a new class of implicit Taylor series methods for 
solving ODE initial value problems. Based on the results presented herein, we conclude 
that the main advantages of implicit series methods are (i) ability to achieve high accuracy 
while maintaining stability, (ii) robustness, (iii) simplicity, and (iv) versatility (ability to 
solve a wide range of problems) . 

However, implicit series methods suffer the disadvantage of being computationally 
intensive. Table XXV summarizes the median number of Newton iterations per step for 
the results presented In the last three sections. The Newton convergence criteria for these 
results was that successive solution iterates differ by less than 10 -14 . 

It is seen that a minimum of four Newton iterations are required per step - two 
iterations for each independent solution calculated. This means that, for a system of 
equations, at least two L-U decompositions must be performed at each step. Consequently, 
an implicit series method must be able to solve a given problem in relatively fewer steps 
to be competitive with other methods. For most nonstiff problems, this implies the need 
for a series with high order accuracy. For stiff problems, this implies the need to vary the 
order of accuracy during the calculation. 

We conclude with the following remarks. First, for nonstiff problems, corrector itera- 
tion is an attractive alternative to Newton iteration. It is easier to implement, requires less 
derivative information, and for small error tolerances is often more computationally effi- 
cient. It works especially well for small, uniform stepsizes. Second, for problems without 
a singular point, Scheme 2T is preferred over Scheme 2. Although both schemes employ 
a degree-two Taylor series, use the same local error estimate (to 0(h 5 )), and generally 
require the same number of marching steps, Scheme 2T can be programmed in real arith- 
metic and is significantly faster as a result. Third, due to the length of this paper, we 
have not made any direct comparisons with other methods. Comparing the efficiency and 
reliability of integration methods has practically become a specialty in itself, and requires 
considerable care and discussion. The intent here is to present enough results to give some 
idea of performance on a range of problems. The next step toward a full comparison with 
other methods is to extend implicit series methods to become a variable order approach. 
This is proposed as a future research topic. Fourth, and lastly, implicit Taylor series meth- 
ods offer a unique combination of high accuracy, stability, robustness, and versatility. The 
author is not aware of another existing method that can, without modification, solve both 
regular ODE’s and stiff systems, together with linear and nonlinear ODE’s with a singular 
point - all with an arbitrarily high order of accuracy. With computer speed and memory 
steadily increasing, algorithm features such as robustness, versatility, and ease of use are 
becoming increasingly important. For this reason, implicit Taylor series methods merit 
consideration as a standard integration tool for solving ordinary differential equations. 
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Problem: 


Scheme 

1 

2 

3 

4 

5 

6a 

6b 

1 

4.00 

5.00 

4.00 

5.00 

4.92 

4.00 

4.00 

2 

4.00 

5.97 

4.12 

5.73 

5.73 

4.00 

4.00 

3 

4.00 

7.72 

4.33 

7.94 

6.90 

4.00 

4.00 

4 

4.00 

7.84 

4.81 

8.56 

7.52 

4.00 

4.00 


9 

10 

11 

12 

13 



1 

4.00 

4.00 

5.95 

6.00 

5.56 



2 

4.00 

4.02 

6.41 

7.01 

5.43 



2T 

4.00 

4.02 

5.82 

6.00 

5.35 



3 

4.02 

4.15 

9.78 

13.65 

6.35 



4 

4.08 

4.18 

17.08 

26.57 

10.56 




Table XXV Median number of Newton iterations per t.imestep. 
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Appendix 


(k) 

Let y n+ 1 := y^-h be the solution of 

h 2 h 3 

yn +1 = Vn + hf(t n+1 ,y n+1 ) - — f'(t n+l: y n+ i) + — /"(f„ +1 ,y n+1 ) 


+ ... + (-l) fc 1 / (fc ^(Wli 2/n+l)- 


(A.l) 


Then for all A; > 2, 


e i+i : = y{tn+ 1) - yi+! 


V ; ’ (k + 1)! yn (lb + 2)! 


(*+l)/? +1) + (* + 2)^/( fc ) 


+ (^7 [ (t+1 f + 2) 4‘ +2) + (t + 1)(* + 3) ^ /<‘ + » 


(A + 2) (A: + 3) df n (fc) 

2 dy /n j 


+ 0(/i fc+4 ). 


(A. 2) 


Proof: 

From equation (A.l), one can expand the solution y^ 1 ^ in a series about the solution 
(k) 

Vn+i f° obtain 


»i* + V ) 


Vn + h /((„+], vh%) + (t/iVi 11 




+ . . . 


+ (-l) fc 


h k+l 

W+Ty- 





- *&) 



< (fc+i) 


(*) \ 2 ^ 2 / (fc) 
^n+ij Q y2 


(An+l?yi+l) “A • 



0/(*) 

dy 


(Li+i» y 


(fc) 
n + 1 


) 


(A.3) 
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from which one gets 



Jk+i) 
Vn + 1 

— Vn + h f(t n+ l,yn+l) - 2/i+l) 


+ . . 

. — . . . 


+ 

(v ik+l) 

yy-n+i 

- v {k) 1 

Vn+ 1 ) 

3/ (fc) x h2 d f u (fc) \ 

1 dy 2 1 (^n+l ? 2/n-f-l) 


+ . . . 

- a>] 

1 , 

+ T 

(ev 1 

- t/ (fc) l 2 

yn+1 j 

S 2 /,. ..(»), 'i 2 0 2 /',. ..(*), 

1 Oyi^n+liVn+l) 21 ( 'n+ljJ/n+l) 


uk+i a2e(k) 


+ ... 


(A.4) 


Now 


Vn + /t/^n+l^i+l) 


e 

2 ! 


f'(tn+l,yn + l) 


+ ... - ... + (- 1 ) l " 1 ir^' ! " ,) (*"+ i-»Si) = i/i+i 


from equation (A.l). Therefore (A.4) becomes 


(A; + l) _ (fc) _ /_ j\fc 
2/n+l ^ri+1 \ 


h fc+1 


/ (fc) ( / 'n+l,yi+l) 


(k + 1)! 

..(fc) x ^ 2 3/', 


+ (2/i+ + i X) - Vn+l) [ / *aj( fn+1 ’ y "+l) ” ^^(Wi.Vn+l) 


+ ... - ... + ( i)» A> +‘ a -^<> - 


(fc + l)! dy 


dy 

(^n+li 2/n+l) 


+ 


1 („(*+!) _ u <*> ^ 

2 \^n+l tfn+1 ^ 


2r , 3 2 / (i (fc) x h 2 d 2 f {k) 

h ^ o (in+l? 2/n+l) 01 /^,2 (^n+lj^rH-l) 


•V 


+ 


+ (_!)* h”^_ay}^, i Jk) 


( k + 1)! dy 2 


2! dy 

(^n+l? Vn+l ) 


+ 


(A. 5) 


(A.6) 


It follows from (A.6) that j/^ 1 ^ - is 0(h k+1 ). Retaining terms through 0(h k+3 ), 

one gets 
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h fe+1 

y': + r - VnU = (-i)‘7rTnT/ (t, ('»+.,!/;5,) 


(fc+1) „(fc) _ 

y, — 

. / (fc+i) (fc) \ f, 9f (fc) v h 2 df‘ (fc) . ^/,3\ 

+ (^n+1 - y„+i J - y-^-^n+l^n+l) + °\ h ) 

so that 


(fc + 1)! 

.2 ar/ 


(fc+1) (fc) 

2/n+l - Vnil 


(-j)* J‘ + Yyi / (t> (<.»+ i,»Si) 


1 - fc^(t„ + ,,!$,) + y + 0(ft 3 ) 


For sufficiently small fc, one then obtains 

x. /i fe+1 


wi+ + i 1} - yi+i = (-i)' 


(A:+l)! /(fc) ( <n+1, ^ t + 1 ) j 1 + 


(A.7) 


(A-8) 


+ 


a/ 

L^y 


^(Wi^Si) 2 - “(Wi,ySi)| |> + o(h fc+4 ). 


(A.9) 


Now each term on the right hand side of (A.9) can be expanded in a series about ( t n , y n ), 
so that 


J fc +D _ Jk) _ 

ifn+l 


y ( n7i = (-i) fe 


h k+ 1 


o 2 fL k) / ( fc) _ \ 

+ 2 at 2 + r n+1 • A v atay 

fy(fc) _ „ ^ ^ 
atay ' \ 2/n+1 dy 2 

For all fc > 2, equation (3.8) shows that 


(fc+1)! 

9 2 /+ 


/n fc) + ^ + (^i+l ~~ 2/n) 


a/, 


(fc) 

n 


dy 




dy 2 


1 + h 


a/n 


/) 2 /■ 

, i 2 u J™ , i 

+ « tttt; + h 


+ fc 2 


r a/»v2 i df n ’ 

ay 


2 ay 


+ 0(h 3 ) 


?/! t +l ~ Vn = hf n + y /' + 0(fc 3 ). 
Substituting (A.ll) into (A. 10) and simplifying, one finally obtains 


yi+V’ - = (- 1 )*+ 


h fc+1 


(fc + 1)! 


/+ + h (/, 


(fc+1) 


+ 


Ofn 

dy 


/+) 


+ 


5! f r(fc+2) . 9 f( fc + 1 ) -+ fW} 

2 V n + dy fn + dy Jn ) 


dy 


dy 

(A.10) 


(A.ll) 


(A.12) 
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Since 


y‘+V - ySi = V «»+.) - - (wt*»+i) - j&V) = «■ 

one obtains from (A. 12) 


(k) _ (fc+i) 

n+1 e n+l 


(A.13) 


,(*+ 1) (fc) 


"n+1 c n+l 


= (-i) 


fc+i 


/i fe+1 

(k + Ty. 


liP + * (/^ +1> + fn k) ) 


dy 


+Y(/r 2) + 2^/r 11 + 


(A- 14) 


Since n and h are both fixed, (A. 14) is a linear difference equation on A;. One can verify by 
direct substitution that the solution to the difference equation is given by (A. 2), plus an 
arbitrary constant which may be taken to be zero. Therefore (A. 2) is valid for all k > 2. 
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Figure l.a Stability region for Euler’s method. 
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Figure 2 Stability region for a degree-3 explicit Taylor series. 
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Stability region for a degree-3 Taylor series with fi = | + i 
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Figure 6 Stability region for a degree-4 fully implicit Taylor series. 
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Figure 7. a Stability region for a degree-4 Taylor series with /i = i + i 
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Stepsize h 

Figure 8.c Reduction of error for Problem 3. Results are shown for the Centered Euler 
Method, Trapezoidal Rule, and Implicit Midpoint Rule. 
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Stepsize h 


Figure 8.e Reduction of error for Problem 5. Results are shown for the Centered Euler 
Method, Trapezoidal Rule, and Implicit Midpoint Rule. 
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Figure 9 Integration path for a degree-2 Taylor series with n — 5 + * g • 
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Figure 12. c 
(J-2 — ^ + i 


Reduction of error for Problem 3 using a degree-3 Taylor series with pi 
A, and with extrapolation. 
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Figure 12. e Reduction of error for Problem 5 using a degree-3 Taylor series with /ij 
/i 2 = | + i and with extrapolation. 
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Figure 13 Stability region for a degree-4 Taylor series with extrapolation 


NASA/TM— 2000-209400 


83 




Max Error 


10 " 
10 ' 
10 ' 
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Figure 14. a Reduction of error for Problem 1 using a degree-4 Taylor series with n\ 

\ + i \ tan jqi H 2 = s + * 5 tan and with extrapolation. 
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Figure 14. b Reduction of error for Problem 2 using a degree-4 Taylor series with Hi 
\ i \ tan yq, H2 = 5 + i § tan fo ’ an< ^ extrapolation. 
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Figure 14. c Reduction of error for Problem 3 using a degree-4 Taylor series with ji\ 

h + i h tan H 2 = \ + * 5 tan fo"- an d with extrapolation. 
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Figure 14. d Reduction of error for Problem 4 using a degree-4 Taylor series with /x i 

\ + i \ tan /x 2 = \ + i | tan and with extrapolation. 
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Figure 14.e Reduction of error for Problem 5 using a degree-4 Taylor series with /x x 
| + i | tan jq, H 2 = | + * I tan W’ anf i with extrapolation. 
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Figure 15. a Stability region for a degree-5 Taylor series with extrapolation. 
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Figure 15. b Stability region for a degree-6 Taylor series with extrapolation. 
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Stepsize h 

Figure 16. a Reduction of error for Problem 1 using a degree-5 Taylor series with fii 
H 2 = h + i h tan and with extrapolation. 
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Figure 17.a Reduction of error for Problems 6.a,b using a degree-1 Taylor series with 
V = \- 


NASA/TM— 2000-209400 


93 




Max Error 



Stepsize h 


Figure 17. b Reduction of error for Problems 6.a,b using a degree-2 Taylor series with 
t l = 2 + * 
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Figure 17.d Reduction of error for Problem 6.b using a degree-3 Taylor series with 
Hi — 4, //2 = \ + i b, and with extrapolation. 
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Figure 17.e Reduction of error for the Problem 6. a using a degree-4 Taylor series with 
/ii = \ + i \ tan yjj, // 2 = \ \ tan and with extrapolation. 
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Figure 19.b Reduction of error for Problems 8.a,b using a degree-2 Taylor series with 
^ = 2 + * : eT- 


NASA/TM— 2000-209400 


104 




tO | H-* 
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Figure 19. c Reduction of error for Problem 8. a using a degree-3 Taylor series with fii 

, = b + i o? an d with extrapolation. 


NASA/TM — 2000-209400 


105 




Max Error 



Stepsize h 


Figure 19. d Reduction of error for Problem 8.b using a degree-3 Taylor series with 
Hi — H2 = i + and with extrapolation. 
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Figure 19. e Reduction of error for Problem 8. a using a degree-4 Taylor series with /ii 

i + i i tan (J - 2 = \ + i § tan and with extrapolation. 
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Figure 19. f Reduction of error for Problem 8.b using a degree-4 Taylor series with /xi 
i 4- i | tan ji- ^2 = ^ + i 5 tan 4^, and with extrapolation. 
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